  ### Results 1.

# - Class Fixed Effects only. This is the main factor I should be
# concerned with.
# - Including SES variables and other variables. I should have up to
# three models and compare then. One with individual characteristics,
# other with indiv. + SES and other controlling for class FE.

ind_controls <- paste("male", "factor(religion_simple)", 
                      "age_grade_distortion", "score_self_esteem", 
                      "score_parents_support", "score_study", sep = "+")

ses_controls <- paste("scores_poverty", "score_neighborhood_quality", 
                      "factor(occ_father_simple_no_job_no_na)",
                      "factor(occ_mother_simple_no_job_no_na)", sep = "+")

class_fe <- paste("factor(class_id)")


regressions_1 <- function(dep_var, race_var, grade_var) {
  # dep_var = "ssi_same_negro_na_2_std"
  # race_var = "factor(negro_na)"
  # grade_var = "grades_score"
  
  # basic regression
  basic_reg <- paste0(dep_var, "~", paste0(race_var,"*", grade_var))
  
  Form_basic <- formula(basic_reg)
  
  # Individual control
  Form_ind <- formula(paste(basic_reg, ind_controls, sep = "+"))
  
  # SES control
  Form_ses <- formula(paste(basic_reg, ses_controls, sep = "+"))

  # Individual + SES control
  Form_ind_ses <- formula(paste(basic_reg, ind_controls, 
                                        ses_controls, sep = "+"))
  
  # Individual controls + Classroom FE
  Form_ind_class_fe <- formula(paste(basic_reg, ind_controls, 
                                          class_fe, sep = "+"))
  
  # SES controls + Classroom FE
  Form_ses_class_fe <- formula(paste(basic_reg, ses_controls, 
                                          class_fe, sep = "+"))
  
  # Ind. + SES controls + Classroom FE
  Form_ind_ses_class_fe <- formula(paste(basic_reg, ind_controls,
                                          ses_controls, class_fe, 
                                          sep = "+"))

  list_specifications <- c(Form_basic, Form_ind, Form_ses, Form_ind_ses,
                           Form_ind_class_fe, Form_ses_class_fe, 
                           Form_ind_ses_class_fe) 
  
  results <- lapply(list_specifications, function(x) (lm(x, data = data, 
                                                         subset = (race < 4))))
  results_r <- lapply(results, robust_se, "HC3")
  results_c <- lapply(results, cluster_se, "class_id")
  
  results_final <- list(results, results_r, results_c)  
  return(results_final)
}

regs_1_s_negro_na <- regressions_1(dep_var = "ssi_same_negro_na_2_std", 
                                   race_var = "negro_na",
                                   grade_var = "grades_score")

regs_1_a_negro_na <- regressions_1(dep_var = "ssi_all_race_2_std", 
                                   race_var = "negro_na", 
                                   grade_var = "grades_score")

regs_1_o_negro_na <- regressions_1(dep_var = "ssi_other_race_negro_na_2_std", 
                                   race_var = "negro_na", 
                                   grade_var = "grades_score")


variables_to_omit <- "(age)|(male)|(religion)|(class_id)|(avg_)|(class)|
                      |(school)|(missing)|(high_school)|(occ_)|
                      |(score_)|(poverty)|(Constant)|(supply)|(n_)|(occ_pai_)"

cov_var_labels_bi <- c("Nonwhite", "Grades", "Nonwhite*Grades")

list_regs_c <- list(regs_1_a_negro_na[[3]][[4]],
                    regs_1_a_negro_na[[3]][[7]],
                    regs_1_s_negro_na[[3]][[4]],
                    regs_1_s_negro_na[[3]][[7]],
                    regs_1_o_negro_na[[3]][[4]],
                    regs_1_o_negro_na[[3]][[7]])

list_regs <- list(regs_1_a_negro_na[[1]][[4]],
                  regs_1_a_negro_na[[1]][[7]],
                  regs_1_s_negro_na[[1]][[4]],
                  regs_1_s_negro_na[[1]][[7]],
                  regs_1_o_negro_na[[1]][[4]],
                  regs_1_o_negro_na[[1]][[7]])

n_obs <- sapply(list_regs, get_obs)

r_adj <- sapply(list_regs, get_adj_r)

testing <- sapply(list_regs, function(x) 
  race_score_linear_test(x, race_inter = "negro_na"))

F_results <- round(c(testing[1,]),3)
P_results <- round(c(testing[2,]),3)

stargazer(list_regs_c, 
          column.separate = c(2,2,2),
          covariate.labels = cov_var_labels_bi,
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other races"),
          omit = variables_to_omit, 
          omit.stat = c("f", "rsq","ser"),
          add.lines = list(c("Linear Hypothesis:", "", "", "", "", "", ""),
                           c("F-statistic:",  F_results),
                           c("P-value", P_results),
                           c("Controls:", "",  "", "",  "",  "",  ""),
                           c("Classroom FE", "No", "Yes", "No", "Yes", "No", "Yes"),
                           c("Observation", n_obs),
                           c("Adjusted R2", r_adj)),
          column.sep.width = "0pt",
          font.size = "tiny", type = "text")
  
  stargazer(list_regs_c, 
            column.separate = c(2,2,2),
            covariate.labels = cov_var_labels_bi,
            dep.var.labels.include = F,
            dep.var.caption = "Social status",
            column.labels = c("All races", "Same race", "Other races"),
            omit = variables_to_omit, 
            omit.stat = c("f", "rsq","ser"),
            add.lines = list(c("Linear Hypothesis:", "", "", "", "", "", ""),
                             c("F-statistic:",  F_results),
                             c("P-value", P_results),
                             c("Controls:", "",  "",  "",  "",  "",  ""),
                             c("Classroom FE", "No", "Yes", "No", "Yes", "No", "Yes"),
                             c("Observation", n_obs),
                             c("Adjusted R2", r_adj)),
            column.sep.width = "-5pt",
            font.size = "tiny", type = "latex", out = "tables/regs_1.tex",
            title = "Academic performance and social status across racial groups",
            label = "regs_1")
  
  

